An actuarial risk model
A classical model for the surplus $U_t$ at time $t$ of an insurance company is the Cramér–Lundberg model (see e.g. Gerber & Shiu (1998)) given by
\[ U_t = U_0 + \gamma t - \sum_{i=1}^{N_t} C_i\]
where $U_0$ is the initial capital, $\gamma$ is a constant premium rate received from the insurees, $C_i$ is a random variable representing the value of the $i$-th claim paid to a given insuree, and $N_t$ is the number of claims up to time $t$. The process $\{N_t\}_t$ is modeled as a Poisson counter, so that the accumulated claims form a compound Poisson process. It is also common to use inhomogeneous Poisson processes and Hawkes self-exciting process, or combinations of such processes for the incidence of the claim, but the classical model uses a homogeneous Poisson counter.
The model above, however, does not take into account the variability of the premium rate received by the company, nor the investiment of the accumulated reserves, among other things. Several diffusion type models have been proposed to account for these and other factors. We will consider a simple model, with a randomly perturbed premium and with variable rentability.
More precisely, we start by rewriting the above expression as the following jump (or impulse) differential equation
\[ \mathrm{d}U_t = \gamma\;\mathrm{d}t - \mathrm{d}C_t,\]
where
\[ C_t = \sum_{i=1}^{N_t} C_i.\]
The addition of an interest rate $r$ leads to
\[ \mathrm{d}U_t = r U_t \mathrm{d}t + \gamma\;\mathrm{d}t - \mathrm{d}C_t.\]
Assuming a premium rate perturbed by a white noise and assuming the interest rate as a process $\{R_t\}_t$, we find
\[ \mathrm{d}U_t = R_t U_t\;\mathrm{d}t + \gamma\;\mathrm{d}t + \varepsilon\;\mathrm{d}W_t - \mathrm{d}C_t,\]
so the equation becomes
\[ \mathrm{d}U_t = (\gamma + R_t U_t)\;\mathrm{d}t + \varepsilon\;\mathrm{d}W_t - \mathrm{d}C_t.\]
Since we can compute exactly the accumulated claims $C_t$, we subtract it from $U_t$ to get rid of the jump term. We also subtract an Ornstein-Uhlenbeck process, in the classical way to transform an SDE into a RODE. So, defining
\[ X_t = U_t - C_t - O_t\]
where $\{O_t\}_t$ is given by
\[ \mathrm{d}O_t = -\nu O_t\;\mathrm{d}t + \varepsilon\;\mathrm{d}W_t,\]
we obtain
\[ \mathrm{d}X_t = (\gamma + R_t U_t)\;\mathrm{d}t + \nu O_t\;\mathrm{d}t = (\gamma + R_t (X_t + C_t + O_t))\;\mathrm{d}t + \nu O_t\;\mathrm{d}t.\]
This leads us to the linear random ordinary differential equation
\[ \frac{\mathrm{d}X_t}{\mathrm{d}t} = R_t X_t + R_t (C_t + O_t) + \nu O_t + \gamma.\]
This equation has the explicit solution
\[ X_t = X_0 e^{\int_0^t R_s\;\mathrm{d}s} + \int_0^t e^{\int_s^t R_\tau\;\mathrm{d}\tau} (R_s (C_s + O_s) + \nu O_s + \gamma)\;\mathrm{d}s.\]
As for the interest rate process $\{R_t\}_t$, there is a vast literature with models for it, see e.g. Chapter 3 of Brigo & Mercurio (2006), in particular Table 3.1. Here, we consider the Dothan model (Section 3.2.2 of the aforementioned reference), which consists simply of a geometric Brownian motion process
\[ \mathrm{d}R_t = \mu R_t \;\mathrm{d}t + \sigma R_t\;\mathrm{d}t,\]
with $R_t = r_0$, where $\mu, \sigma, r_0$ are positive constants. This has an explicit solution
\[ R_t = r_0 e^{(\mu - \sigma^2/2)t + \sigma W_t},\]
so that the above equation for $\{X_t\}_t$ is a genuine random ODE.
Once the solution of $\{X_t\}_t$ is obtained, we find an explicit formula for the surplus $X_t = U_t - C_t - O_t$, namely
\[ U_t = C_t + O_t + X_0 e^{\int_0^t R_s\;\mathrm{d}s} + \int_0^t e^{\int_s^t R_\tau\;\mathrm{d}\tau} (R_s (C_s + O_s) + \nu O_s + \gamma)\;\mathrm{d}s,\]
with $\{R_t\}_t$ as above.
Numerical simulations
Setting up the problem
First we load the necessary packages:
using Plots
using Random
using LinearAlgebra
using Distributions
using RODEConvergenceThen we define the random seed:
rng = Xoshiro(123)Random.Xoshiro(0xfefa8d41b8f5dca5, 0xf80cc98e147960c1, 0x20e2ccc17662fc1d, 0xea7a7dcb2e787c01)The evolution law:
function f(t, x, y)
o = y[1]
r = y[2]
c = y[3]
ν = 5.0
γ = 1.0
dx = r * (x + c + o) + ν * o + γ
return dx
endf (generic function with 1 method)The time interval:
t0, tf = 0.0, 3.0(0.0, 3.0)The law for the initial condition:
x0 = 1.0
x0law = Dirac(x0)Distributions.Dirac{Float64}(value=1.0)The Ornstein-Uhlenbeck, geometric Brownian motion, and compound Poisson processes for the noise term:
O0 = 0.0
Oν = 5.0
Oε = 0.8
R0 = 0.2
Rμ = 0.02
Rσ = 0.4
Cmax = 0.2
Cλ = 8.0
Claw = Uniform(0.0, Cmax)
noise = ProductProcess(
OrnsteinUhlenbeckProcess(t0, tf, O0, Oν, Oε),
GeometricBrownianMotionProcess(t0, tf, R0, Rμ, Rσ),
CompoundPoissonProcess(t0, tf, Cλ, Claw)
)ProductProcess{Float64, Tuple{OrnsteinUhlenbeckProcess{Float64}, GeometricBrownianMotionProcess{Float64}, CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}}}((OrnsteinUhlenbeckProcess{Float64}(0.0, 3.0, 0.0, 5.0, 0.8), GeometricBrownianMotionProcess{Float64}(0.0, 3.0, 0.2, 0.02, 0.4), CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}(0.0, 3.0, 8.0, Distributions.Uniform{Float64}(a=0.0, b=0.2))))The resolutions for the target and approximating solutions, as well as the number of simulations for the Monte-Carlo estimate of the strong error:
ntgt = 2^18
ns = 2 .^ (6:9)
nsample = ns[[1, 2, 3, 4]]
m = 400400And add some information about the simulation:
info = (
equation = "a risk model",
noise = "coupled Ornstein-Uhlenbeck, geometric Brownian motion, and compound Poisson processes",
ic = "\$X_0 = $x0\$"
)(equation = "a risk model", noise = "coupled Ornstein-Uhlenbeck, geometric Brownian motion, and compound Poisson processes", ic = "\$X_0 = 1.0\$")We define the target solution as the Euler approximation, which is to be computed with the target number ntgt of mesh points, and which is also the one we want to estimate the rate of convergence, in the coarser meshes defined by ns.
target = RandomEuler()
method = RandomEuler()RandomEuler{Float64, Distributions.Univariate}(Float64[])Order of convergence
With all the parameters set up, we build the ConvergenceSuite:
suite = ConvergenceSuite(t0, tf, x0law, f, noise, target, method, ntgt, ns, m)ConvergenceSuite{Float64, Distributions.Dirac{Float64}, ProductProcess{Float64, Tuple{OrnsteinUhlenbeckProcess{Float64}, GeometricBrownianMotionProcess{Float64}, CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}}}, typeof(Main.var"Main".f), 2, 1, RandomEuler{Float64, Distributions.Univariate}, RandomEuler{Float64, Distributions.Univariate}}(0.0, 3.0, Distributions.Dirac{Float64}(value=1.0), Main.var"Main".f, ProductProcess{Float64, Tuple{OrnsteinUhlenbeckProcess{Float64}, GeometricBrownianMotionProcess{Float64}, CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}}}((OrnsteinUhlenbeckProcess{Float64}(0.0, 3.0, 0.0, 5.0, 0.8), GeometricBrownianMotionProcess{Float64}(0.0, 3.0, 0.2, 0.02, 0.4), CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}(0.0, 3.0, 8.0, Distributions.Uniform{Float64}(a=0.0, b=0.2)))), RandomEuler{Float64, Distributions.Univariate}(Float64[]), RandomEuler{Float64, Distributions.Univariate}(Float64[]), 262144, [64, 128, 256, 512], 400, [1, 1, 1, 1], [9.3953855e-316 NaN NaN; 1.177204523e-315 NaN NaN; … ; NaN 4.0e-323 0.0; 0.0 NaN 0.0], [0.0, 1.35966498e-315, 1.477686316e-315, 0.0, 0.0, 0.0, 1.35966498e-315, 1.477686513e-315, 0.0, 0.0 … 0.0, 0.0, NaN, 1.328929276e-315, 1.488047347e-315, 0.0, 0.0, NaN, 1.328929276e-315, 1.488047544e-315], [0.0, 9.32698786e-316, 0.0, 0.0, 0.0, 9.32699063e-316, 9.32699063e-316, 0.0, 9.32699063e-316, 0.0 … 9.32718707e-316, 0.0, 9.32718707e-316, 0.0, 0.0, 0.0, 9.32718984e-316, 9.32718984e-316, 0.0, 9.32718984e-316])Then we are ready to compute the errors via solve:
@time result = solve(rng, suite) 3.774541 seconds (299.78 k allocations: 19.306 MiB, 12.12% compilation time)The computed strong error for each resolution in ns is stored in result.errors, and a raw LaTeX table can be displayed for inclusion in the article:
table = generate_error_table(result, info) \begin{tabular}[htb]{|r|l|l|l|}
\hline N & dt & error & std err \\
\hline \hline
64 & 0.0469 & 0.223 & 0.0167 \\
128 & 0.0234 & 0.115 & 0.00878 \\
256 & 0.0117 & 0.0579 & 0.00405 \\
512 & 0.00586 & 0.0286 & 0.00209 \\
\hline
\end{tabular}
\bigskip
\caption{Mesh points (N), time steps (dt), strong error (error), and standard error (std err) of the Euler method for a risk model for each mesh resolution $N$, with initial condition $X_0 = 1.0$ and coupled Ornstein-Uhlenbeck, geometric Brownian motion, and compound Poisson processes, on the time interval $I = [0.0, 3.0]$, based on $M = 400$ sample paths for each fixed time step, with the target solution calculated with $262144$ points. The order of strong convergence is estimated to be $p = 0.987$, with the 95\% confidence interval $[0.9019, 1.0727]$.}The calculated order of convergence is given by result.p:
println("Order of convergence `C Δtᵖ` with p = $(round(result.p, sigdigits=2)) and 95% confidence interval ($(round(result.pmin, sigdigits=3)), $(round(result.pmax, sigdigits=3)))")Order of convergence `C Δtᵖ` with p = 0.99 and 95% confidence interval (0.902, 1.07)Plots
We plot the rate of convergence with the help of a plot recipe for ConvergenceResult:
plt_result = plot(result)For the sake of illustration of the behavior of the system, we visualize a sample solution
plt_sols = plot(suite, ns=nothing, label="\$X_t\$", linecolor=1)We also illustrate the convergence to a sample solution
plt_suite = plot(suite)We can also visualize the noises associated with this sample solution:
plt_noises = plot(suite, xshow=false, yshow=true, label=["\$O_t\$" "\$R_t\$" "\$C_t\$"], linecolor=[1 2 3])The actual surplus is $U_t = X_t - O_t - C_t$, so we may visualize a sample solution of the surplus by subtracting these two noises from the solution of the above RODE.
plt_surplus = plot(range(t0, tf, length=ntgt+1), suite.xt .- suite.yt[:, 1] .- suite.yt[:, 3], xaxis="\$t\$", yaxis="\$u\$", label="\$U_t\$", linecolor=1)This page was generated using Literate.jl.